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An  Introduction  to  Singular  Perturbations  "Sy— — — 

The  books  of  Wasow  (1965),  Vasil 'eva  and  Butuzov  (1973),  and  O'Malley 
(1974)  all  discuss  the  asymptotic  solution  of  singularly  perturbed  initial 


value  problems  for  vector  systems  of  the  form 


Gypc  t 


a  =  f(a,3,t,c),  a(0)  prescribed 
e3  =  g(a,3,t,e),  3(0)  prescribed 


/lt!  3  4- 


on  a  finite  or  semi-infinite  interval  where  t  ^  0.  Here  c  is  a  small  positive 
parameter  and  an  asymptotic  solution  valid  as  c  -*•  0  is  sought.  Assuming 
a  little  smoothness  and  assuming  that  the  eigenvalues  of  the  Jacobian  matrix 
g  are  strictly  stable  (i.e.,  have  strictly  negative  real  parts)  everywhere 


implies  that  the  asymptotic  solution  has  the  form 


I 


a(t)  =  M(t)  +  0(c) 

B(t)  =  N(t)  +  n(t/e)  +  0(e) 


t 

\3 


where  M  and  N  satisfy  the  "reduced"  system 


M  =  f(M,N,t,0),  M(0)  =  a(0) 
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and  n(x)  is  an  asymptotically  stable  solution  of  the  (parameterless) 
"boundary  layer"  problem 


(4) 


=  g(M(0),  N(0)  +  n»0,0),  t  >  0 
n(0)  =  3(0)  -  N(0). 


The  0(e)  terms  represent  small  residuals  bounded  in  magnitude  (on  fixed 
finite  t  intervals)  by  Ce,  for  some  bounded  constant  C,  provided  e  is 
sufficiently  small.  Hoppensteadt  (1966)  further  showed  that  such  resultsTjfc? 
continue  to  hold  for  all  t  ^  0  provided  M  and  N  decay  exponentially  as 
t 

It  is  essential  to  realize  the  substantial  order  reduction  achieved 


3d 
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y 
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through  the  approximations  (2).  The  implicit  function  theorem  guarantees' on> 


a  locally  unique  solution 


(5) 


N(t)  =  $(M(t),t) 


®iet. 
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of  the  nonlinear  system  g(M,N,t,0)  =  0.  (In  practice,  numerical  or 
symbolic  solution  methods  would  be  used  to  obtain  <j>,  especially  when  3 
has  many  components  and  the  order  reduction  is  substantial.)  Knowing  <j>, 
M(t)  follows  from  solving  the  initial  value  problem 


(6) 


M  =  F(M.t)  s  f(M,4>(M,t),t,0),  M(0)  =  a(0). 


With  our  hypotheses  (or,  e.g.,  the  weaker  ones  of  Howes  (1979)),  a  unique 

r 

solution  h(t)  of  (4)  will  decay  to  zero  as  the  "stretched  variable"  t/e  -*• 

Thus,  the  n-term  allows  nonuniform  convergence  (as  e  -*■  0)  of  the  3  variable 
in  a  narrow  "boundary  layer"  region  near  t  =  0.  Since  n  becomes  negligible  as 
e  ■*  0  for  each  fixed  t  >  0,  the  limiting  solution  away  from  t  =  0  is  uniquely 
determined  by  the  solution  M(t)  of  the  reduced  problem  (6)  and  by  N(t)  =  4>(M(t),t). 
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The  reduction  of  order  is  completely  analogous  to  the  pseudo-steady 
state  hypothesis  of  chemical  kinetics.  In  the  celebrated  enzyme  kinetics 
theory  of  Michael  is  and  Menton,  an  e  can  be  introduced  (as  in  (1))  as  the 
ratio  of  the  initial  enzyme  to  the  initial  substrate  concentrations  (cf. 
Heinekin  et  al.  (1967)  and  Murray  (1977)).  Solving  the  reduced  problem, 
instead  of  the  original  problem  with  its  fast  initial  transient  (repre¬ 
sented  by  n(t/e)  in  (2))  is  most  advantageous  for  purposes  of  numerical 
integration  (cf.  Aiken  and  Lapidus  (1974)  and  Miranker  (1980)).  The  full 
system  will  be  stiff  for  small  values  of  e,  and  stability  requirements  for 
Euler's  method  and  similar  explicit  schemes  will  force  one  to  use  a  small 
stepsize  for  t  >  0,  even  though  the  solution  is  smooth  there.  As  Shampine 
(1980)  and  Dahlquist  et  al . (1 980)  point  out,  implicit  schemes  for  the  full 
system  may  not  so  restrict  the  stepsize  and  may  actually  be  solving  systems 
closely  linked  to  the  reduced  problem  (3). 

The  preceding  hypotheses  are  restrictive,  and  it  should  be  realized 
that  related  (though  less  complete)  results  are  available  in  other  (mere 
complicated)  contexts  (cf.,  e.g.,  Bavinck  and  Grasman  (1969),  Hoppensteadt 
and  Miranker  (1976),  O'Malley  and  Flaherty  (1980),  Bobisud  and  Christenson 
(1980),  and  O'Malley  (1980)).  If,  for  example,  gy  is  (everywhere)  singular, 
we  cannot  expect  to  uniquely  determine  the  limiting  solution  from  the  re¬ 
duced  problem  alone.  If  g^  has  purely  imaginary  eigenvalues,  we  can  expect 
rapidly  oscillating  solutions.  If  has  eigenvalues  with  both  positive  and 
negative  real  parts,  the  boundary  layer  system  is  conditionally  stable  and  we 
can  expect  nonuniform  convergence  (boundary  layer  behavior)  for  corresponding 
two-point  problems  at  both  endpoints.  Considerable  complexity  can  be  expected 
when  (where)  the  eigenvalue  structure  of  g^  changes  (both  bifurcation  theory 
and  turning  point  theory  may  be  called  for). 


Real  problems  naturally  simultaneously  involve  many  small  parameters. 


Detailed  theories  are  available  for  problems  with  several  interrelated  para¬ 
meters  (cf.  O'Malley  (1969))  and,  even,  for  t-dependent  parameters  (cf. 

Gingold  (1980)).  The  one  parameter  e  in  (1)  should  be  interpreted  as  an 
aggregate  parameter  used  to  distinguish  the  "fast"  variables  8  from  the  "slow" 
variables  a.  More  refined  models  could  use  several  small  parameters  to 
distinguish  different  groupings  of  fast  variables.  Indeed,  a  realistic 
modeler  often  introduces  a  hierarchy  of  models,  increasing  in  dimensionality, 
allowing  increasingly  accurate  models  of  the  same  system.  In  practice,  re¬ 
searchers  will  make  different  decisions  regarding  both  the  overall  number  of 
components  employed  and  the  appropriate  number  of  variables  for  reduced- 
order  models. 

The  structure  of  the  singular  perturbation  problem  (1)  is  very  special, 
because  the  components  a  are  clearly  separated  (via  the  parameter  e)  from 
the  6  components  which  are  anticipated  to  be  more  rapidly  varying,  at  least 
in  narrow  boundary  or  interior  transition  regions.  In  applied  problems, 
practitioners  often  know  which  variables  are  "fast."  There  is,  for  example, 
a  natural  dichotomy  when  one  deals  with  dynamics  involving  both  electrical 
and,  relatively  sluggish,  mechanical  components.  Analogous  obvious  splittings 
commonly  occur  in  chemical  contexts  as  well.  Often,  however,  variables  are 
not  completely  decoupled  into  fast  and  slow  sets  and  a  preliminary  decoupling 
procedure  must  be  used  before  singular  perturbation  concepts  become  helpful 
(i.e.,  before  a  meaningful  small  parameter  e  can  be  introduced.)  Numerical 
decoupling  methods  are  currently  being  developed  by  Enright  and  Kamel  and  by 
Dahlquist  and  Sbderlind  (preliminary  work  is  reported  in  Enright  and  Kamel 
(1979)  and  Soderlind  (1979)).  We  also  note  the  emphasis  of  Dahlquist  et  al. 
(1980)  on  scaling  techniques  for  nonlinear  differential  systems  in  dimension¬ 
less  variables  (cf.  Lin  and  Segel  (1974)  as  well).  Our  experience  relates 
primarily  to  models  for  aircraft  engines,  structural  dynamics,  and  power 
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distribution  and  generation  (cf.  O'Malley  and  Anderson  (1980),  Anderson 
and  Hallauer  (1980),  Kokotovic  et  al.(1980),  and  Girijashankar  et  al. 
(1980)).  We  expect,  however,  that  the  underlying  strategy  will  be  appli¬ 
cable  in  chemical  systems  and  a  wide  spectrum  of  other  fields. 

Linear  Systems  without  Explicit  Parameters 

We  shall  consider  large-scale  linear  systems  of  the  form 

(7)  x  =  A(t)x  +  u(t) 

where  u(t)  is  a  forcing  function,  not  a  feedback  control.  Such  a  system 
could,  of  course,  be  obtained  from  linearizing  a  nonlinear  system  about  a 
nominal  trajectory.  The  initial  behavior  of  any  solution  will  depend  on 
the  eigenvalues  of  A(0).  Let  us  suppose  that 

(HI).  n1  of  the  eigenvalues  of  A(0)  are  small 

in  magnitude  compared  to  the  remaining 
n2  =  n-n^  eigenvalues. 

(A  finer  subdivision  may  be  appropriate  for  a  more  careful  model.)  If 
such  a  slow/fast  split  in  solution  modes  for  the  homogeneous  system  could 
be  maintained,  we  might  hope  (as  for  the  singularly  perturbed  system  (1)) 
that  the  approximate  dynamics  of  the  full  n  th  order  system  (7)  (away  from 
transient  regions)  could  be  given  by  a  lower  th  order  system. 

We  will  first  try  to  transform  system  (7)  to  a  decoupled  form  by  use 
of  a  transformation 

(8)  y  =  Ht)x  =  T2(t)  Tj(t)x  =  r2(t)z 
where  and  T2  are  the  block-triangular  matrices 
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specified  by  the  decoupling  matrices  L ( t )  and  K(t).  Since  these  matrices 
are  easily  inverted,  we  immediately  have  the  inverse 


-L  I  +  LK  . 

"2 

If  we  subdivide  A  =  (A..)  after  its  first  n,  rows  and  columns,  the  inter- 

I J  I 

mediate  variable  z  =  T^x  will  satisfy 

z  =  B(t)z  +  7^  (t)u 

with  the  compatibly  partitioned  B  *  (B^)  being  upper  block  triangular 
provided  L(t)  satisfies  the  matrix  Riccati  equation 

L  =  A22  L  -  LAn  +  LA12  L  - 

A21 


(ID 


Specifically,  we  note  that 


*  ! 
i 


(12)  B.j ^  =  A-j i  -  ®12  =  ^12*  ^21  ~  ^22  =  ^22  "*"  ^12* 


We  also  note  that  transformations  such  as  (9)  occur  in  many  classical  con¬ 
texts,  including  reduction-of-order  (cf.  Hartman  (1964))  and  asymptotic 
theory  (cf.  Wasow  (1965)).  Here,  the  matrix  L  is  generally  nonsquare  and 
the  Riccati  equation  has  time-varying  coefficients,  in  contrast  to  the 
symmetric  and  time-invariant  case  most  familiar  in  control  and  estimation 
theory  (cf.  Van  Dooren  (1980)).  In  the  general  context  needed  here,  con¬ 
ditions  guaranteeing  existence  of  the  Riccati  solution  do  not  seem  to  be 
known.  Most  available  theory  is  contained  in  Reid  (1972),  Medanic  (1979), 
and  Bart  et  al.  (1979).  Blowup  of  solutions  can,  of  course,  be  detected 
in  numerical  solutions  and  "integration  to  blowup"  is  actually  used  to 
calculate  eigenlengths  (cf.  Nelson  and  Elder  (1977)).  It  would  be  helpful 
to  know  which  solution  of  the  Riccati  equation  one  should  use  (i.e.,  which 
initial  value  to  impose  for  L(0) )  to  prolong  the  interval  of  existence  or 
to  optimize  the  stability  of  a  numerical  solution  technique. 

Using  the  second  transforming  matrix  T^it)  and  splitting  the  variable  y 
after  its  first  n^  rows  provides  the  decoupled  system 


(13) 


>22 


r  \ 

*i 


+  T(t)u 


provided  the  matrix  K  satisfies  the  linear  system 


K  8  B„  K  -  KB,,  -  A 


'11 


22  "12. 


(14) 


(A  block  diagonal ization  is  actually  unnecessary  to  decouple  modes  of 
the  unforced  system  since  already  achieved  a  block  triangularization, 
but  we  proceed  to  diagonalize  because  obtaining  K  is  relatively  easy.) 

Altogether,  then,  we  have  shown  that  solving  the  original  vector  or 
matrix  system  (7)  (thereby  finding  a  fundamental  matrix  for  the  homo¬ 
geneous  problem)  is  equivalent  to  solving  the  four  systems 

L  =  A22L  -  LA^  +  LA12  L  -  A21 
K  =  (An  -  A12  L)  K  -  K(A22'+  LA12)  -  A12 

(16) 

*i  =  (An  '  Ai2  l)  yi  +  (Ini  +  KL  K)u 

y2  =  (A22  +  LA12^  y2  +  !n2^u 

We  hope  to  convince  the  reader,  that  with  the  right  hypotheses,  we  can 
determine  solutions  to  these  systems  easier  than  solving  the  original 
system  directly. 

Note  that  if  L(0)  =  0,  7^(0)  will  be  a  similarity  transformation, 
and  B(0)  and  C(0)  will  have  the  same  eigenvalues  as  A(0).  We  would  have 
succeeded  in  separating  the  initially  fast  (and  slow)  modes  of  the  trans¬ 
formed  homogeneous  system  if  we  could  guarantee  that  the  n-|  smallest 
eigenvalues  of  B(0)  (identified  by  hypothesis  (HI))  are  those  of  8^(0) 

(the  n2  largest  eigenvalues  will  then  automatically  be  those  of  B22(0)). 

The  algebraic  Riccati  equation  L (0)  =  0  for  L(0)  will  have  many  solutions, 
but  that  solution  L(0)  which  makes  the  eigenvalues  of  8^(0)  =  A^(0)>A12(0)  L(0) 

coincide  with  the  n^  small  eigenvalues  of  A(0)  is  unique.  This  result  is 
well-known  for  the  symmetric  algebraic  Riccati  equation  (cf.  Martensson 


(1971))  and  it  was  proved  in  Anderson  (1979)  in  our  nonsquare  con¬ 
text.  Anderson  also  showed  that  the  appropriate  L(0)  could  be 
readily  obtained  by  iterating.  For  computational  purposes,  one 
defines 

(17)  *i+1  »  (A 22(0)  +  £1  A12(0))_1  (Zi  An(0)  +  A12(0)) 

for  each  i  >  0  and  lets  L(0)  =  lim  £...  This  scheme  turns  out  to 

*j  CO 

be  robust  with  respect  to  initial  iterates  &  and  it  converges  like 
a  geometric  series  with  ratio  u  proportional  to  the  quotient  of  the 
largest  small  to  the  smallest  large  eigenvalue  of  A(0).  An  alter- 

-1  ~M1~ 

native  representation  is  L(0)  *  -M,  M,  where  spans  the  n,  di- 

c  1  M 

m2 

L  J 

mensional  (generalized)  eigenspace  of  the  small  eigenvalues  of  A(0). 

(A  less  convenient  computational  approach  to  finding  L(0)  would  then 
be  to  approximate  a  basis  for  this  eigenspace  and  form  the  indicated 
quotient.)  This  formula  makes  it  clear  that  a  rearrangement  of  the 
components  of  x  may  be  necessary  to  assure  the  invertibility  of  .  Any 
ill-conditioning  of  can  also  cause  difficulties  (cf.  Watkins  (1980)). 

We  shall  assume  that  our  system  is  "two-time  scale"  after  the  trans¬ 
formation,  viz 

(H2)  the  solution  L(t)  of  the  initial  value  problem  for 
the  matrix  Riccati  equation  remains  bounded  through¬ 
out  a  fixed  interval  0  £t£  T  and  the  eigenvalues  of 
the  matrix  B^(t)  =  -  A^2  L  remain  small  in  mag¬ 

nitude  compared  to  those  of  B22(t)  =  A22  +  L  A^2 


throughout  the  interval . 


I  -  r"  n  ■  i  rlg*l 


Under  this  hypothesis,  the  free  (i.e. ,  unforced)  response  of  the  first 
n^  components  of  the  decoupled  system  for  y  will  be  much  slower  than 
that  of  the  last  n2  components.  {For  two  point  problems,  T  will  represent 
the  right  endpoint,  while  it  should  be  a  typical  length  scale  for  initial 
value  problems.)  In  practice,  it  may  be  advisable  to  reinitialize  L(t) 
from  time-to-time,  just  as  the  Scott  and  Watts’  (1977)  boundary  value  code 
reorthogonal izes  to  maintain  stability. 

Specializing  to  initial-value  problems  for  (7),  let  us  also  impose  a 
"fast-mode  stability"  hypothesis,  viz 


(H3)  the  eigenvalues  of  B22(t)  all  have  large  strict!' 
negative  real  parts  throughout  0  1  t  £  T. 


In  using  this  assumption,  we  are  close  to  succumbing  to  a  popular  heresy: 
that  eigenvalue  stability  implies  stability.  This  is  incorrect  for  time- 
varying  systems,  as  Hoppensteadt  (1966)  points  out  by  counterexample.  As 
our  initial  statement  shows,  however,  it  is  true  for  singularly  perturbed 
systems  with  appropriately  smooth  coefficients.  Kreiss  (1978,  1979) 
presents  a  counterexample,  but  not  in  the  special  structure  generally 
assumed.  We  could  carry  out  an  explicit  and  rigorous  analysis  involving 
the  two  small  parameters  implicitly  introduced  by  assumptions  (H2)  and  (H3), 
but  this  would  require  definite  assumptions  about  parameter  dependence 
that  we  wish  to  avoid  being  explicit  about. 

We  note  that  the  differential  system  (14)  for  K  is  opposite  in  stability 
to  the  linearized  equation 


l  =  +  £  B22 


(18) 


T 
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for  L.  Since  (H2)  implies  that  the  initial  value  problem  for  L(t)  is 
well-behaved,  (backward)  numerical  integration  of  a  terminal  value 
problem  for  K  is  suggested.  Taking  K(T)  =  o  uniquely  determines  K(T) 
since  the  resulting  system  for  K(T)  is  a  linear  one  whose  coefficient 
matrices  ^ (T)  and  B 22 (T)  have  no  eigenvalues  in  common  (cf.,  e.g.. 

Bellman  (1970)). 

Because  (H-3)  implies  that  B ^ (t)  has  eigenvalues  with  large  nega¬ 
tive  real  parts,  the  system  (14)  for  K(t)  can  be  rewritten  in  the  trad¬ 
itional  singular  perturbation  form  with  a  small- parameter  multiplying 
the  derivative.  Indeed,  the  usual  theory  implies  that  the  limiting 
solution  for  t  <  T  will  satisfy  a  reduced  system  K(t)  =  0  (i.e.  a  pseudo- 
steady  state  approximation  will  apply.).  Because  we  have  K(T)  =  0,  we'd 
expect  no  discontinuity  at  T.  K(t)  can,  then,  be  readily  calculated  as 
the  limit  of  the  iteration 

(19)  Kjt,  -  -B-2\  (A,2  -  B„  Kj) 

(The  derivative  term  can  also  be  retained  as  a  small  perturbation  in  the 
right-hand-side.)  Thus,  the  differential  system  (14)  for  K  has  been 
essentially  replaced  by  an  easier-to-solve  linear  algebraic  equation 
K(t)  =  0.  In  conducting  numerical  experiments,  we've  found  that  the 
pseudo-steady  state  algebraic  Riccati  equation  L(t)  *  0  is  quite  a  good 
approximation  for  finding  L(t)  for  t  >  0,  but  we  are  not  willing  to  generally 
recommend  that  drastic  simplification.  We  note  that  the  linearization  (18) 
(about  a  solution  l  of  our  Riccati  system)  is  singularly  perturbed  under 
our  hypotheses,  but  the  original  system  (11)  is  not  of  singular  perturbation 
form.  SSderlind  and  Dahlquist  (personal  conversation)  are,  however,  ex¬ 
periencing  success  in  using  piecewise  constant  solutions  to  algebraic 
Riccati  equations  in  similar  contexts. 


Once  L  and  K  are  (approximately)  determined,  there  remain  linear 
systems  for  y-|  and  y^.  The  "fast"  system 

(20)  y2  =  B22(t)y2  +  [L  I  ]  u(t) 

will  be  amenable  to  a  pseudo-steady  state  approximation  y2  =  0  away  from 
t  =  0  provided  the  resulting  approximate  y2s  is  slowly-varying.  Thus, 


we'll  finally  ask 


(H4)  The  approximation 

*2s(t>  -  -B  -2\  (L  ln)  ii. 
is  slowly- varying 

(i.e.  y2s  is  small)  throughout  0  <  t  <_  T.  Even  though  ’s  lar9e,  note 
that  this  might  be  difficult  to  achieve  if  u(t)  were  rapidly  varying.  If 
the  initial  behavior  of  y2  is  important,  one  needs  to  add  a  "fast"  initial 
layer  corrector  t0  the  "slow-mode"  solution  y2s-  The  corrector  will 
be  a  decaying  solution  to  =  B22(t)y2f..  Because  of  our  hypothesis  (H3), 
this  system  would  only  need  to  be  integrated  -  with  a  small  stepsize  - 
over  a  short  initial  interval.  This  correction  will  be  initially  important 
because  y2$(0)  ?  y2(°)  will  generally  force  an  initial  nonuniform  convergence 
in  y2  components. 

The  remaining  "slow"  system  (16)  for  y-j  will  be  integrated  as  an  initial 
value  problem  using  the  transformed  initial  condition 

y,(0)  =  [I  +  K(0)L(0)  K( 0) ]  x(0).  Having  found  (or  approximated)  L,  K,  y, , 
i  n  j  * 

and  y2,  the  inverse  transformation 

T  'n.  -K(l>  1  h(t) 


■L(t)  In^+  L(t)  K(t)  y2(t) 


provides  (an  approximation  for)  the  original  variable  x(t).  Altogether, 
for  t  >  0,  our  approximations  allow  us  to  obtain  an  approximate  funda¬ 
mental  matrix  for  x  by  integrating  the  systems  for  L  and  y-j  and  using 
algebraic  systems  for  K  and  y2- 

In  practice,  we've  solved  initial  value  problems  for  vectors  of 
dimensionsfour  to  sixteen  (cf.  O'Malley  and  Anderson  (1980))  with  two 
to  five  corresponding  slow  modes  providing  order  reduction.  Away  from 
the  initial  point,  good  approximations  were  soon  achieved  by  the  low 
order  models  even  though  the  parameters  measuring  time-scale  separation 
and  fast-mode  stability  were  only  moderately  small.  For  such  low  order 
models,  we  have  the  advantage  of  being  able  to  compare  results  to  accurate 
solutions  obtained  numerically  with  small  stepsizes.  More  testing  is 
needed  for  higher  dimensional  models,  where  it  is  essential  to  use  one's 
physical  insight  to  a  maximum  to  make  initial  decisions  on  which  components 
are  expected  to  be  predominantly  slow  or  fast.  More  effort  should  also  be 
spent  on  boundary  value  problems.  There  the  advantages  of  our  asymptotic 
techniques  will  be  even  greater  than  for  initial  value  problems. 

Use  of  some  of  these  concepts  for  nonlinear  problems  has  been  reported  by 
Kreiss  (1979),  Girijashankar  et  al.  (1980),  Kokotovic  et  al.  (1979),  and 
Watkins  (1980). 
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